clear all;
%% Code explanation
% 180902
% 190710 updated to obtain Tc Bc2 with 50%Rn, 1%Rn, added offset to R1 R2
% 200121 DZ corrected the resistance factors, export to spreadsheet data
%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter
L_W1=43/66; % 3Sn_13PbTe
L_W2=95/61; % 5Sn_15PbTe
ratio=0.5;

%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
root=('D:\Documents2\Rawdata\Sn\180801 MPI mK\20180820_inplane\'); %

filenames={'31_0uW_roots';'29_20uW_roots';'28_40uW_roots';'27_60uW_roots';'26_80uW_roots';'25_100uW_roots';...
    '24_140uW_roots';'23_200uW_roots';'22_300uW_roots';'21_400uW_roots';'20_600uW_roots';'20_800uW_roots';...
    '19_180uW';'19_300uW';'18_500uW';'18_800uW';'18_1100uW';'17_1500uW';'15_2000uW';'15_2500uW';'15_3200uW';...
    '15_4000uW';'13_4000uW';'13_4500uW';'13_5000uW';'13_5500uW';'12_6000uW';'11_7000uW'};



 Rnorm=[1.36/L_W1/10 2.4/L_W2/10]'; % normal state sheet resistances of the two samples
NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,2);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT+1,2);

Tmax1=0.6818;
Tmax2=1.219;


for kk=1:NT
    file1=[root,'201808',filenames{kk},'.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f %f','commentStyle','#');
        data1=data1';  
        fclose(fid);
    
     
          B=round(data1{1}*1000)/1000;
          R1=(data1{2}/L_W1+0.002)/10;
           R2=(data1{3}/L_W2+0.005)/10;
          T=T_RuOx(data1{5}*100);
          Tmc_err(kk,1)=3*std(T);
          
         B0=(min(B):0.02:max(B))';
        NB=1;
        for jj=2:size(B,1)
            if B(jj)~=B(jj-1)
                NB=NB+1;
            end
        end
        Bset=zeros(NB,1);
        jj2=1;
        Bset(1)=B(1);
        for jj=2:size(B,1)
            if B(jj)~=B(jj-1)
               jj2=jj2+1;
               Bset(jj2)=B(jj);
            end
        end
         RT1=zeros(length(Bset),2); % T and R1
         RT2=zeros(length(Bset),2); % T and R2
    
         
         for ii=1:length(Bset)
            indexB=find(B==Bset(ii));
            [RT1(ii,2),j1]=min(R1(indexB));
            RT1(ii,1)=T(indexB(j1));
            [RT2(ii,2),j2]=min(R2(indexB));
            RT2(ii,1)=T(indexB(j2));
         end
          RRT1=interp1(Bset,RT1(:,2),B0,'spline');
          TT1=interp1(Bset,RT1(:,1),B0,'linear');
         
          RRT2=interp1(Bset,RT2(:,2),B0,'spline');
          TT2=interp1(Bset,RT2(:,1),B0,'linear');
          
          subplot(2,2,1),plot(Bset,RT1(:,2),'s','Color',colorset(kk,:)','MarkerSize',3);
          hold on
          subplot(2,2,3),plot(B0,RRT1,'-','Color',colorset(kk,:)','LineWidth',1);
          hold on
           subplot(2,2,2),plot(Bset,RT2(:,2),'s','Color',colorset(kk,:)','MarkerSize',3);
          hold on
          subplot(2,2,4),plot(B0,RRT2,'-','Color',colorset(kk,:)','LineWidth',1);
          hold on

          
          indexHc=find(RRT1<=ratio*Rnorm(1));
          if isempty(indexHc)~=1
              Hc2(kk,1)=B0(indexHc(length(indexHc)));
              Tmc(kk,1)=TT1(indexHc(length(indexHc)),1);
          end
          indexHc=find(RRT2<=ratio*Rnorm(2));
          if isempty(indexHc)~=1
              Hc2(kk,2)=B0(indexHc(length(indexHc)));
              Tmc(kk,2)=TT2(indexHc(length(indexHc)),1);
          end
            title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[Bset RT2(:,2) RT2(:,1) mean(RT2(:,1))*ones(length(Bset),1)];
       
         xlswrite('FigS5B.xlsx',AA,kk);
       xlswrite('FigS5B.xlsx',title,kk);
%     
end
Bplot=[0:0.1:4]';
subplot(2,2,1),plot(Bplot,ratio*Rnorm(1)*ones(length(Bplot)),'-r');
hold on
subplot(2,2,2),plot(Bplot,ratio*Rnorm(2)*ones(length(Bplot)),'-r');
hold on

figure
plot([Tmc(:,1)' Tmax1]',Hc2(:,1),'dr',[Tmc(:,2)' Tmax2]',Hc2(:,2),'sb');
hold on
% 



%% END




